
clearvars
starttime = datetime('now');

% ////////////////////////////////////////////////////////
% SET UP FILES AND PARAMETERS
% ////////////////////////////////////////////////////////

main_base = " ";

%bname = "strainup15_comm";
nrun = 5;
fnames = strings(1,nrun);
names = ["n2_AP", "n4_AP", "n5_AP", "n6_AP", "n9_AP"];
for i=1:nrun
    
    %single level corona reanalysis
    fnames(i) = strcat(main_base, names(i));
    %if i>3
    %    fnames(i) = strcat(main_base, "whole embryos/",names(i));
    %else
    %    fnames(i) = strcat(main_base, "free explants/", names(i));
    %end
    
end



for q=1:nrun %1:nrun
    
    
    javaaddpath 'C:\Program Files\MATLAB\R2019b\java\mij.jar'
    javaaddpath 'C:\Program Files\MATLAB\R2019b\java\ij.jar'
    MIJ.start('C:\Users\Lab\Fiji.app')
    basedir = fnames(q);
    
    %basedir = "C:/Users/Lab/Box Sync/Backups/Macros Sommer/seg_test/ext_amp/seven";
    
    % write base file for sim to read
    fid=fopen('C:\Users\Lab\Downloads\basefile.txt','wt');
    fprintf(fid, '%s', basedir');
    fclose(fid);
    
    
    % ////////////////////////////////////////////////////////
    % MACRO RUNNING
    % ////////////////////////////////////////////////////////
    
    cd("C:\Users\Lab\OneDrive - University of Pittsburgh\_BoxMigration\Backups\Macros Sommer\seg_test")
    %run only t1 track
    % MIJ.run("t1track live")
    
    % %{
    % //comment block start for just running sim and segmentation, no analysis
    %these need to be installed to the local FIJI as plugins
    MIJ.run("cellstrain new box")
    MIJ.run("corona live basef onelv CP")
    MIJ.run("t1track live CP")
    
    % ////////////////////////////////////////////////////////
    % ANALYSIS
    % ////////////////////////////////////////////////////////
    
    % ********* cell and corona strain analysis **************
    
    % Import Cellstrain
    % Script for importing data from the following text file:
    %
    %    filename: C:\Users\Lab\Box Sync\Backups\Macros Sommer\seg_test\CellEllipse.txt
    %
    % Auto-generated by MATLAB on 05-Nov-2021 13:48:08
    opts = delimitedTextImportOptions("NumVariables", 4);
    
    % Specify range and delimiter
    opts.DataLines = [2, Inf];
    opts.Delimiter = "\t";
    
    % Specify column names and types
    opts.VariableNames = ["SliceNumA", "CellID", "CentroidX", "CentroidY", "Angle", "major", "minor", "Area", "Width", "Height", "Perim"];
    opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double"];
    
    % Specify file level properties
    opts.ExtraColumnsRule = "ignore";
    opts.EmptyLineRule = "read";
    
    % Import the data
    CellEllipse = readtable(strcat(basedir, "/CellEllipse.txt"), opts);
    
    
    % Clear temporary variables
    clear opts
    
    % Import CoronaData
    % Script for importing data from the following text file:
    %
    %    filename: C:\Users\Lab\Box Sync\Backups\Macros Sommer\seg_test\Coronas\CoronaData.txt
    %
    % Auto-generated by MATLAB on 05-Nov-2021 16:16:34
    
    % Setup the Import Options and import the data
    opts = delimitedTextImportOptions("NumVariables", 12);
    
    % Specify range and delimiter
    opts.DataLines = [2, Inf];
    opts.Delimiter = "\t";
    
    % Specify column names and types
    opts.VariableNames = ["CellID", "SliceNumA", "CoronaCells", "CentroidX", "CentroidY", "Angle", "major", "minor", "Width", "Height", "Area", "sided"];
    opts.VariableTypes = ["double", "double", "string", "double", "double", "double", "double", "double", "double", "double", "double", "double"];
    
    % Specify file level properties
    opts.ExtraColumnsRule = "ignore";
    opts.EmptyLineRule = "read";
    
    % Specify variable properties
    opts = setvaropts(opts, "CoronaCells", "WhitespaceRule", "preserve");
    opts = setvaropts(opts, "CoronaCells", "EmptyFieldRule", "auto");
    
    % Import the data
    CoronaData = readtable(strcat(basedir, "\Coronas\CoronaData.txt"), opts);
    
    
    % Clear temporary variables
    clear opts
    
    % Strain analysis
    
    % set p = IDs of zapped cells to speed up if extrusion
    
    % Setup the Import Options and import the data
    opts = delimitedTextImportOptions("NumVariables", 1);
    
    % Specify range and delimiter
    opts.DataLines = [2, Inf];
    opts.Delimiter = ",";
    
    % Specify column names and types
    opts.VariableNames = "Value";
    opts.VariableTypes = "double";
    
    % Specify file level properties
    opts.ExtraColumnsRule = "ignore";
    opts.EmptyLineRule = "read";
    
    
    % Clear temporary variables
    clear opts
    %
    cd("C:/Users/Lab/OneDrive - University of Pittsburgh/_BoxMigration/Backups/Macros Sommer/seg_test/");
    
    if (isa(CellEllipse,'table'))
        CellEllipse = table2struct(CellEllipse);
    end
    
    if (isa(CoronaData,'table'))
        CoronaData = table2struct(CoronaData);
    end
    
    
    %     cells2 = unique([CoronaData.CellID]);
    %
    %
    %
    %       a=[CellEllipse([CellEllipse.SliceNumA]==1).CellID];
    %       b=[CoronaData([CoronaData.SliceNumA]==1).CellID];
    %       c=setdiff(b,a);
    %      [bn,ibn,ib]=intersect(cells,c);
    %       cells(ibn)=[];
    
    
    uniqueSliceNums1 = unique([CellEllipse.SliceNumA]);
    uniqueCellIDs1 = unique([CellEllipse.CellID]);
    
    uniqueSliceNums2 = unique([CoronaData.SliceNumA]);
    uniqueCellIDs2 = unique([CoronaData.CellID]);
    
    % Find CellIDs present in all slice numbers for each table
    idsInAllSlices1 = [];
    for i = 1:length(uniqueCellIDs1)
        currentID = uniqueCellIDs1(i);
        q = [CellEllipse.CellID];
        r = [q==currentID];
        idSliceNums = [CellEllipse(r).SliceNumA];
        %idSliceNums = CellEllipse.SliceNumA(CellEllipse.CellID == currentID);
        if isequal(sort(unique(idSliceNums)), sort(uniqueSliceNums1))
            idsInAllSlices1 = [idsInAllSlices1; currentID];
        end
    end
    
    idsInAllSlices2 = [];
    for i = 1:length(uniqueCellIDs2)
        currentID = uniqueCellIDs2(i);
        q = [CoronaData.CellID];
        r = [q==currentID];
        idSliceNums = [CoronaData(r).SliceNumA];
        %idSliceNums = CoronaData.SliceNumA(CoronaData.CellID == currentID);
        if isequal(sort(unique(idSliceNums)), sort(uniqueSliceNums2))
            idsInAllSlices2 = [idsInAllSlices2; currentID];
        end
    end
    
    % Find common CellIDs present in all slice numbers in both tables
    commonIDsInAllSlices = intersect(idsInAllSlices1, idsInAllSlices2);
    cells = commonIDsInAllSlices';
    
    % uniqueTimeFrames = unique([CoronaData.SliceNumA]);
    % uniqueIDs = unique([CoronaData.CellID]);
    
    
    % Initialize a logical array to track presence of each ID
    % idsPresentInAllFrames = {};
    %
    % % Loop through each ID
    % for i = 1:length(uniqueIDs)
    %     currentID = uniqueIDs(i);
    %     % Check if the current ID is present in each time frame
    %     for j = 1:length(uniqueTimeFrames)
    %         q = [CoronaData.CellID];
    %         r = [q==currentID];
    %         idTimeFrames = [CoronaData(r).SliceNumA];
    %         % Logical indexing to check presence
    %          if isequal(sort(unique(idTimeFrames)), sort(uniqueTimeFrames))
    %             idsPresentInAllFrames{end+1} = currentID; % Append to cell array
    %          end
    %     end
    % end
    %
    % % Extract IDs present in all time frames
    % idsPresentInAllFrames = uniqueIDs(idInAllTimeFrames);
    % cells = idsPresentInAllFrames;
    
    
    % endpoint strain calculation
    [CellStrain_end,~] = CellStrainCalc_endpoint(CellEllipse, 20, cells);
    [CoronaStrain_end,~] = CellStrainCalc_endpoint(CoronaData,20, cells);
    
    [CellStrain,~] = CellStrainCalc(CellEllipse,cells, 20);
    [CoronaStrain,~] = CellStrainCalc(CoronaData,cells,20);
    
    
    
    
    % Get data out
    
    
    if (isa(CellEllipse,'table'))
        CellEllipse = table2struct(CellEllipse);
    end
    
    if (isa(CoronaData,'table'))
        CoronaData = table2struct(CoronaData);
    end
    % Allocating Variables up front
    %cells = unique([CoronaData(:).CellID]); %finds all the unique corona IDs
    % ^cells, usually
    
    time =20; % update based on number of frames
    
    
    cd("C:/Users/Lab/OneDrive - University of Pittsburgh/_BoxMigration/Backups/Macros Sommer/seg_test/");
    
    
    G = zeros(time,length(cells));
    GG = zeros(time,length(cells));
    K = zeros(time,length(cells));
    KK = zeros(time,length(cells));
    M = zeros(time,length(cells));
    MM = zeros(time,length(cells));
    g = zeros(time, 1);
    gg = zeros(time, 1);
    k = zeros(time, 1);
    kk = zeros(time, 1);
    m = zeros(time, 1);
    mm = zeros(time, 1);
    
    
    for jj=1:length(cells)
        % g is corona strain in the xx direction.
        % gg is cell strain in the xx direction
        
        % k is corona strain Eyy
        % kk is cell strain in the yy direction
        %cells(jj)
        
        % m is Area strain of corona
        % mm is Area strain of cell
        
        for i=2:time
            %   i
            %collect Exx of that cell's corona
            g(i)=[CoronaStrain_end([CoronaStrain_end.CellID]==cells(jj) & [CoronaStrain_end.time]==i).Exx];
            
            
            %Doing the same as above for the YY direction:
            k(i)=[CoronaStrain_end([CoronaStrain_end.CellID]==cells(jj) & [CoronaStrain_end.time]==i).Eyy];
            
            
            gg(i)=[CellStrain_end([CellStrain_end.CellID]==cells(jj)& [CellStrain_end.time]==i).Exx];
            
            kk(i)=[CellStrain_end([CellStrain_end.CellID]==cells(jj)& [CellStrain_end.time]==i).Eyy];
            
            %AC is area strain stuff
            mm(i)= [CellStrain_end([CellStrain_end.CellID]==cells(jj)& [CellStrain_end.time]==i).EA];
            m(i) = [CoronaStrain_end([CoronaStrain_end.CellID]==cells(jj)& [CoronaStrain_end.time]==i).EA];
            
        end
        
        G(:,jj) = g;
        GG(:,jj) = gg;
        K(:,jj) = k;
        KK(:,jj) = kk;
        M(:,jj) = m;
        MM(:,jj) = mm;
        
    end
    
    S = zeros(time,length(cells));
    SS = zeros(time,length(cells));
    V = zeros(time,length(cells));
    VV = zeros(time,length(cells));
    U = zeros(time,length(cells));
    UU = zeros(time,length(cells));
    s = zeros(time, 1);
    ss = zeros(time, 1);
    v = zeros(time, 1);
    vv = zeros(time, 1);
    u = zeros(time, 1);
    uu = zeros(time, 1);
    
    for jj=1:length(cells)
        % g is corona strain in the xx direction.
        % gg is cell strain in the xx direction
        
        % k is corona strain Eyy
        % kk is cell strain in the yy direction
        %cells(jj)
        
        % m is Area strain of corona
        % mm is Area strain of cell
        
        for i=2:time
            %   i
            %collect Exx of that cell's corona
            s(i)=[CoronaStrain([CoronaStrain.CellID]==cells(jj) & [CoronaStrain.time]==i).Exx];
            
            
            %Doing the same as above for the YY direction:
            v(i)=[CoronaStrain([CoronaStrain.CellID]==cells(jj) & [CoronaStrain.time]==i).Eyy];
            
            
            ss(i)=[CellStrain([CellStrain.CellID]==cells(jj)& [CellStrain.time]==i).Exx];
            
            vv(i)=[CellStrain([CellStrain.CellID]==cells(jj)& [CellStrain.time]==i).Eyy];
            
            %AC is area strain stuff
            uu(i)= [CellStrain([CellStrain.CellID]==cells(jj)& [CellStrain.time]==i).EA];
            u(i) = [CoronaStrain([CoronaStrain.CellID]==cells(jj)& [CoronaStrain.time]==i).EA];
            
        end
        
        S(:,jj) = s;
        SS(:,jj) = ss;
        V(:,jj) = v;
        VV(:,jj) = vv;
        U(:,jj) = u;
        UU(:,jj) = uu;
        
    end
    
    
    % ********* Markov analysis **************
    
    concordx = zeros(length(cells),time);
    concordy = zeros(length(cells),time);
    
    
    % Finding concordance for Exx:
    % Iterate over cells and time to populate the concordx matrix
    for jj = 1:length(cells)
        
        for i = 1:time
            %  Find Exx for the Corona for the cell that the loop is on, for
            %  the timestep that the loop is on
            Cor = [CoronaStrain([CoronaStrain.CellID]==cells(jj) & [CoronaStrain.time]==i).Exx];
            
            % Find Exx for the cell & time that the loop is on
            Cell = [CellStrain([CellStrain.CellID]==cells(jj)& [CellStrain.time]==i).Exx];
            
            % Check if the sign of Exx for the corona is the same as the cell
            % in question at time t
            if sign(Cor) == sign(Cell)
                % If the sign is the same, call this concordance and assign
                % value 1
                concordx(jj,i) = 1;
            else
                % Let 0 signify discordance for that cell and time if the sign
                % does not match
                concordx(jj,i) = 0;
            end
        end
        
    end
    
    
    % find the concordance for Eyy with new matrix name
    
    
    % Iterate over cells and time to populate the concord1 matrix
    for jj = 1:length(cells)
        
        for i = 1:time
            %  Find Exx for the Corona for the cell that the loop is on, for
            %  the timestep that the loop is on
            Cor = [CoronaStrain([CoronaStrain.CellID]==cells(jj) & [CoronaStrain.time]==i).Eyy];
            
            % Find Exx for the cell & time that the loop is on
            Cell = [CellStrain([CellStrain.CellID]==cells(jj)& [CellStrain.time]==i).Eyy];
            
            % Check if the sign of Exx for the corona is the same as the cell
            % in question at time t
            if sign(Cor) == sign(Cell)
                % If the sign is the same, call this concordance and assign
                % value 1
                concordy(jj,i) = 1;
            else
                % Let 0 signify discordance for that cell and time if the sign
                % does not match
                concordy(jj,i) = 0;
            end
        end
        
    end
    
    P1x =  zeros(length(cells),time);
    P2x =  zeros(length(cells),time);
    P3x =  zeros(length(cells),time);
    P4x =  zeros(length(cells),time);
    
    for j = 1:length(cells)
        y_obs = concordx(j,:)';
        y_obs(y_obs == 1) = 2;
        y_obs(y_obs < 1) = 1;
        L = length(y_obs);
        T = L;
        
        P_MC = zeros(2,2);
        
        for t=1:L-1
            P_MC(y_obs(t),y_obs(t+1))= P_MC(y_obs(t),y_obs(t+1))+1;
            
            P1x(j,t) =  P_MC(1,1)/sum(P_MC(1,:));
            P2x(j,t) =  P_MC(1,2)/sum(P_MC(1,:));
            P3x(j,t) =  P_MC(2,1)/sum(P_MC(2,:));
            P4x(j,t) =  P_MC(2,2)/sum(P_MC(2,:));
        end
        
    end
    
    P1y =  zeros(length(cells),time);
    P2y =  zeros(length(cells),time);
    P3y =  zeros(length(cells),time);
    P4y =  zeros(length(cells),time);
    
    for j = 1:length(cells)
        y_obs = concordy(j,:)';
        y_obs(y_obs == 1) = 2;
        y_obs(y_obs < 1) = 1;
        L = length(y_obs);
        T = L;
        
        P_MC = zeros(2,2);
        
        for t=1:L-1
            P_MC(y_obs(t),y_obs(t+1))= P_MC(y_obs(t),y_obs(t+1))+1;
            
            P1y(j,t) =  P_MC(1,1)/sum(P_MC(1,:));
            P2y(j,t) =  P_MC(1,2)/sum(P_MC(1,:));
            P3y(j,t) =  P_MC(2,1)/sum(P_MC(2,:));
            P4y(j,t) =  P_MC(2,2)/sum(P_MC(2,:));
        end
        
    end
    
    % ********* signed strain analysis **************
    
    sstrainx = zeros(length(cells),time);
    
    for jj = 1:length(cells)
        
        for i = 2:time
            %  Find Exx for the Corona for the cell that the loop is on, for
            %  the timestep that the loop is on
            Cor = [CoronaStrain([CoronaStrain_end.CellID]==cells(jj) & [CoronaStrain.time]==i).Exx];
            
            % Find Exx for the cell & time that the loop is on
            Cell = [CellStrain([CellStrain.CellID]==cells(jj)& [CellStrain.time]==i).Exx];
            
            % Check if the sign of Exx for the corona is the same as the cell
            % in question at time t
            if sign(Cor) == sign(Cell)
                % If the sign is the same, call this concordance and assign
                % value 1
                sstrainx(jj,i) = abs(Cell - Cor);
                
            else
                % Let 0 signify discordance for that cell and time if the sign
                % does not match
                sstrainx(jj,i) = -(abs(Cell-Cor)) ;
            end
        end
        
    end
    
    
    sstrainy = zeros(length(cells),time);
    
    for jj = 1:length(cells)
        
        for i = 2:time
            %  Find Exx for the Corona for the cell that the loop is on, for
            %  the timestep that the loop is on
            Cor = [CoronaStrain([CoronaStrain.CellID]==cells(jj) & [CoronaStrain.time]==i).Exx];
            
            % Find Exx for the cell & time that the loop is on
            Cell = [CellStrain([CellStrain.CellID]==cells(jj)& [CellStrain.time]==i).Exx];
            
            % Check if the sign of Exx for the corona is the same as the cell
            % in question at time t
            if sign(Cor) == sign(Cell)
                % If the sign is the same, call this concordance and assign
                % value 1
                sstrainy(jj,i) = abs(Cell - Cor);
                
            else
                % Let 0 signify discordance for that cell and time if the sign
                % does not match
                sstrainy(jj,i) = - ( abs(Cell-Cor) );
            end
        end
        
    end
    
    
    % ********* Shape index analyis and area things  **************
    
    %CellEllipse = table2struct(CellEllipse);
    
    si = zeros(length(cells), time);
    areas = zeros(length(cells), time);
    ars = zeros(length(cells), time);
    angs = zeros(length(cells), time);
    
    for j=1:length(cells)
        for i=1:time
            perim = CellEllipse([CellEllipse.SliceNumA]==i & [CellEllipse.CellID]==cells(j)).Perim;
            area = CellEllipse([CellEllipse.SliceNumA]==i & [CellEllipse.CellID]==cells(j)).Area;
            si(j,i) = perim/sqrt(area);
        end
    end
    
    for j=1:length(cells)
        for i=1:time
            
            area = CellEllipse([CellEllipse.SliceNumA]==i & [CellEllipse.CellID]==cells(j)).Area;
            areas(j,i) = area;
            ang = CellEllipse([CellEllipse.SliceNumA]==i & [CellEllipse.CellID]==cells(j)).Angle;
            angs(j,i) = ang;
            
            major = CellEllipse([CellEllipse.SliceNumA]==i & [CellEllipse.CellID]==cells(j)).major;
            minor = CellEllipse([CellEllipse.SliceNumA]==i & [CellEllipse.CellID]==cells(j)).minor;
            
            ars(j,i) = major/minor;
            
        end
    end
    
    
    % ********* MSD analysis **************
    
    % Use x and y centroids from CellEllipse
    
    
    % Calculate MSDs
    
    % remove first zero, create a new matrix with each cell and time from
    % long array
    
    xCens = zeros(length(cells), time);
    yCens = zeros(length(cells), time);
    
    for i=1:length(cells)
        for j=1:time
            xCens(i,j)=CellEllipse([CellEllipse.SliceNumA]==j & [CellEllipse.CellID]==cells(i)).CentroidX;
            yCens(i,j)=CellEllipse([CellEllipse.SliceNumA]==j & [CellEllipse.CellID]==cells(i)).CentroidY;
        end
    end
    
    
    MSD = zeros(length(cells), time);
    
    for i = 1:length(cells) % over cells
        for j = 1:time  % over time
            sums(i,j) =  ( ( xCens(i,j) - xCens(i,1) )^2 + (  yCens(i,j) -  yCens(i,1) )^2 );
        end
    end
    
    cs = cumsum(sums,2);
    asum = cumsum(areas,2);
    
    for i = 1:length(cells)
        for j = 1:time
            MSD(i,j) = ( (1/j) .* cs(i,j) )/(asum(i,j) );
        end
    end
    
    
    
    % ********* T1 analysis **************
    
    % import data
    % Setup the Import Options and import the data
    opts = delimitedTextImportOptions("NumVariables", 6);
    
    % Specify range and delimiter
    opts.DataLines = [2, Inf];
    opts.Delimiter = ["\t", ","];
    
    % Specify column names and types
    opts.VariableNames = ["t1time", "t1cells", "VarName3", "VarName4", "VarName5", "VarName6"];
    opts.VariableTypes = ["double", "double", "double", "double", "double", "double"];
    
    % Specify file level properties
    opts.ExtraColumnsRule = "ignore";
    opts.EmptyLineRule = "read";
    
    % Import the data
    t1s = readtable(strcat(basedir,"\t1s.txt"), opts);
    
    % Convert to output type
    t1s = table2array(t1s);
    
    % Clear temporary variables
    clear opts
    
    %remove duplicate entries, get rid of NaNs
    t1s(isnan(t1s))=0;
    t1s = unique(t1s,'rows');
    
    %create table of binary events by cell and time
    t1 = zeros(length(cells), time);
    
    %t1 = zeros(length(t1s), time); for testing
    for i=1:size(t1s,1)
        for j=2:size(t1s,2)
            if ismember(t1s(i,j), cells)
                cell_ind = find(cells==t1s(i,j));
                %cell_ind = t1s(i,j); for testing
                t1(cell_ind, t1s(i,1)) = 1;
                
            end
        end
    end
    
    % ********* SIDEDNESS ANALYSIS **************
    
    %doing initial sidedness of cell by finding n cells in one level corona
    %using that value for all values so the info carries over in data
    
    sided = zeros(length(cells), time);
    %CellEllipse([CellEllipse.SliceNumA]==j & [CellEllipse.CellID]==cells(i)).CentroidX;
    for j=1:length(cells)
        for i=1:time
            sided(j,i) = CoronaData([CoronaData.SliceNumA]==i & [CoronaData.CellID]==cells(j)).sided;
            
        end
    end
    
    
    % ////////////////////////////////////////////////////////
    % GENERATE 3D MATRIX WITH VARIABLES OF INTEREST
    % ////////////////////////////////////////////////////////
    
    % first index is data, second is cell, third is time
    
    data = [];
    
    for i=1:length(cells)
        for j=1:time
            
            data(1,i,j) = cells(i); % cell ID
            
            % populate cell and corona strains
            data(2,i,j) = GG(j,i); % cell exx
            data(3,i,j) = KK(j,i); % cell eyy
            data(4,i,j) = G(j,i); % cor exx
            data(5,i,j) = K(j,i); % cor eyy
            data(6,i,j) = MM(j,i); % cell e_area
            data(7,i,j) = M(j,i); % cor e_area
            
            % populate state transition probabilities
            data(8,i,j) = P1x(i,j); % pdd x
            data(9,i,j) = P2x(i,j); % pdc x
            data(10,i,j) = P3x(i,j); % pcd x
            data(11,i,j) = P4x(i,j); % pcc x
            data(12,i,j) = P1y(i,j); % pdd y
            data(13,i,j) = P2y(i,j); % pdc y
            data(14,i,j) = P3y(i,j); % pcd y
            data(15,i,j) = P4y(i,j); % pcc y
            
            % populate signed difference in cell and corona strain
            data(16,i,j) = sstrainx(i,j); % cell-cor x
            data(17,i,j) = sstrainy(i,j); % cell-cor y
            % add area after com
            
            %populate shape index
            data(18,i,j) = si(i,j);
            
            %populate MSD
            data(19,i,j) = MSD(i,j);
            
            %populate T1s
            data(20,i,j) = t1(i,j);
            
            %populate parameters
            %data(21,i,j) = ps_now(1); %viscosity
            %data(22,i,j) = ps_now(2); %packing
            %data(22,i,j) = ps_now(3); %border/corner
            %data(23,i,j) = ps_now(4); %repulsion strength
            %data(24,i,j) = ps_now(5); %attraction strength
            %data(25,i,j) = ps_now(6); %aniso attrac (angle range)
            %data(26,i,j) = ps_now(7); %polarity module
            %data(27,i,j) = ps_now(8); %size rectangle analysis
            %data(28,i,j) = ps_now(9); %attraction rectangle only
            % data(29,i,j) = ps_now(10); %room to grow
            
            %STEP STRAINS
            data(29,i,j) = SS(j,i); % cell exx
            data(30,i,j) = VV(j,i); % cell eyy
            data(31,i,j) = S(j,i); % cor exx
            data(32,i,j) = V(j,i); % cor eyy
            data(33,i,j) = UU(j,i); % cell e_area
            data(34,i,j) = U(j,i); % cor e_area
            
            data(35,i,j) = xCens(i,j);% xc
            data(36, i, j) = yCens(i,j);% yc
            data(37,i,j) = areas(i,j); % area cell
            data(38, i ,j) = ars(i,j);% aspect ratio major/minor
            data(39,i,j)= angs(i,j);% angle
            data(40,i,j) = sided(i,j); %sidedness
            
            
            
        end
    end
    
    
    cd(basedir);
    save('wspace.mat');
    
    %}
    
    %comment block end for just running sim and segmentation, no analysis
    
    
    % ////////////////////////////////////////////////////////
    % POST ANALYSIS END MATTER
    % ////////////////////////////////////////////////////////
    MIJ.run("Clear Results");
    MIJ.closeAllWindows
    MIJ.exit
    javarmpath 'C:\Program Files\MATLAB\R2019b\java\mij.jar'
    javarmpath 'C:\Program Files\MATLAB\R2019b\java\ij.jar'
    
    
    clearvars -except fnames starttime notes
    
end
endtime = datetime('now');
endtime-starttime